{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "# Network Causality : Model-free Inference\n",
    "\n",
    "**Content creator**: Kai Chen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "# Tutorial Objectives\n",
    "\n",
    "In this notebook, we'll start to explore some of the model-free inference methods. More importantly, we gonna to compare performances of those methods. \n",
    "\n",
    "1. Introduction to time-delayed mutual informaiton.\n",
    "1. Causal inference in linear and nonlinear systems.\n",
    "1. Causal inference for experimental data.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "#### References:\n",
    "- D. Zhou, Y. Xiao, Y. Zhang, Z. Xu and D. Cai, “Granger causality network reconstruction of conductance-based integrate-and-fire neuronal systems”, PLoS ONE, 9 (2), e87636, 2014. ([PDF](https://ins.sjtu.edu.cn/people/zdz/publication_papers/Granger_Causality_Reconstruction.pdf))\n",
    "- S. Li, Y. Xiao, D. Zhou and D. Cai, “Causal inference in nonlinear systems: Granger causality versus time-delayed mutual information”, Phys. Rev. E, 97 (5), 052216, 2018. ([PDF](https://ins.sjtu.edu.cn/people/zdz/publication_papers/Granger_causality_versus_time_delayed_mutual_information.pdf))\n",
    "- Matlab-based fast GC-estimator: [GC_clean](https://github.com/bewantbe/GC_clean)\n",
    "- Python-based mutual-information estimator: [mutual-information](https://github.com/NeoNeuron/mutual_information)\n",
    "- SciPy $\\chi^2$ statistics [documents](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "---\n",
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.stats import chi2\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "#@title Figure settings\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "nma_style = {\n",
    "    'figure.figsize' : (8, 6),\n",
    "    'figure.autolayout' : True,\n",
    "    'font.size' : 15,\n",
    "    'xtick.labelsize' : 'small',\n",
    "    'ytick.labelsize' : 'small',\n",
    "    'legend.fontsize' : 'small',\n",
    "    'axes.spines.top' : False,\n",
    "    'axes.spines.right' : False,\n",
    "    'xtick.major.size' : 5,\n",
    "    'ytick.major.size' : 5,\n",
    "}\n",
    "for key, value in nma_style.items():\n",
    "    plt.rcParams[key] = value\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @title Helper functions\n",
    "\n",
    "def create_structure_array(x:np.ndarray, order:int)->np.ndarray:\n",
    "    '''\n",
    "    Prepare structure array for regression analysis.\n",
    "\n",
    "    Args:\n",
    "    x         : original time series\n",
    "    order     : regression order\n",
    "\n",
    "    Return:\n",
    "    x_array   : structure array with shape (len(x)-order) by (order).\n",
    "\n",
    "    '''\n",
    "    N = len(x) - order\n",
    "    x_array = np.zeros((N, order))\n",
    "    for i in range(order):\n",
    "        x_array[:, i] = x[-i-1-N:-i-1]\n",
    "    return x_array\n",
    "\n",
    "def auto_reg(x, order)->np.ndarray:\n",
    "    '''\n",
    "    Auto regression analysis of time series.\n",
    "\n",
    "    Args:\n",
    "    x         : original time series\n",
    "    order     : regression order\n",
    "\n",
    "    Return:\n",
    "    res       : residual vector\n",
    "\n",
    "    '''\n",
    "    reg_array = create_structure_array(x, order)\n",
    "    coef = np.linalg.lstsq(reg_array, x[order:], rcond=None)[0]\n",
    "    res = x[order:] - reg_array @ coef\n",
    "    return res\n",
    "\n",
    "def joint_reg(x, y, order)->np.ndarray:\n",
    "    '''\n",
    "    Joint regression analysis of time series.\n",
    "\n",
    "    Args:\n",
    "    x         : original time series 1\n",
    "    y         : original time series 2\n",
    "    order     : regression order\n",
    "\n",
    "    Return:\n",
    "    res       : residual vector\n",
    "\n",
    "    '''\n",
    "    reg_array_x = create_structure_array(x, order)\n",
    "    reg_array_y = create_structure_array(y, order)\n",
    "    reg_array = np.hstack((reg_array_x, reg_array_y))\n",
    "    coef = np.linalg.lstsq(reg_array, x[order:], rcond=None)[0]\n",
    "    res = x[order:] - reg_array @ coef\n",
    "    return res\n",
    "\n",
    "def GC(x, y, order):\n",
    "    '''\n",
    "    Granger Causality from y to x\n",
    "\n",
    "    Args:\n",
    "    x         : original time series (dest)\n",
    "    y         : original time series (source)\n",
    "    order     : regression order\n",
    "\n",
    "    Return:\n",
    "    GC_value  : residual vector\n",
    "\n",
    "    '''\n",
    "    res_auto = auto_reg(x, order)\n",
    "    res_joint = joint_reg(x, y, order)\n",
    "    GC_value = 2.*np.log(res_auto.std()/res_joint.std())\n",
    "    return GC_value\n",
    "\n",
    "def GC_SI(p, order, length):\n",
    "    '''\n",
    "    Significant level of GC value.\n",
    "\n",
    "    Args\n",
    "    p       : p-value\n",
    "    order   : parameter of chi^2 distribution\n",
    "    length  : length of data.\n",
    "\n",
    "    Return:\n",
    "    significant level of null hypothesis (GC\n",
    "        between two independent time seies)\n",
    "\n",
    "    '''\n",
    "    return chi2.ppf(1-p, order)/length"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "# Section 1: Time-delayed Mutual Information (TDMI)\n",
    "\n",
    "TDMI is an information-theoretic approach for detecting causal interactions. In general, the quantity of mutual information characterizes the common information shared between two signals. Given two stationary signals X and Y with their time traces {$x_t$} and {$y_t$}, respectively, the mutual information between them is defined as\n",
    "$$\n",
    "I(X, Y)=\\sum_{x_{t}} \\sum_{y_t} p\\left(x_{t}, y_{t}\\right) \\log \\frac{p\\left(x_{t}, y_{t}\\right)}{p\\left(x_{t}\\right) p\\left(y_{t}\\right)}\n",
    "$$\n",
    "where $ p(xt ,yt )$ is the joint probability distribution of $X = x_t$ and $Y = y_t$ , $p(x_t)$ and $p(y_t)$ are their marginal probability distributions. In particular, $I (X,Y) = 0$ is equivalent to $p(x_t ,y_t) = p(x_t)p(y_t)$, which indicates that signals $X$ and $Y$ are independent if they do not share information.\n",
    "\n",
    "Mutual information is symmetric, i.e., $I (X,Y ) = I (Y,X)$. Therefore, it cannot be applied directly to infer the direction of interactions between two signals. To overcome this limitation, one can introduce a time-lag parameter $\\tau$ to capture the delay of information transfer between the two signals. TDMI as a function of time-lag $\\tau$ is defined as\n",
    "\n",
    "$$\n",
    "I(X, Y, \\tau)=\\sum_{x_{t}} \\sum_{y_{t-\\tau}} p\\left(x_{t}, y_{t-\\tau}\\right) \\log \\frac{p\\left(x_{t}, y_{t-\\tau}\\right)}{p\\left(x_{t}\\right) p\\left(y_{t-\\tau}\\right)}\n",
    "$$\n",
    "\n",
    "where $p(x_t ,y_{t −\\tau} )$ is the joint probability distribution of $X = x_t$ and $Y = y_{t −\\tau}$ . A nonzero amplitude of the mutual information as a function of $\\tau$ indicates the existence of interactions between two signals; the sign of the time-lag $\\tau$ where $I(X,Y,\\tau)$ reaches its peak magnitude is used to infer the information flow direction that can be further interpreted as the causal direction of interaction. A negative $\\tau$ indicates that $X$ shares a maximum amount of information with the future of $Y$ , thus $X$ drives $Y$. A positive $\\tau$ indicates that $X$ shares a maximum amount of information with the past of $Y$, thus $X$ is driven by $Y$ ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# Add you code to estimate mutual information\n",
    "def mi(x, y, bins = None):\n",
    "    \"\"\"mutual information for 0-1 binary time series\n",
    "    :param x: first series\n",
    "    :type x: int of ndarray\n",
    "    :param y: second series\n",
    "    :type y: int of ndarray\n",
    "    :return: mutual information\n",
    "\n",
    "    \"\"\"\n",
    "    N = len(x)\n",
    "    if bins is None:\n",
    "        bins = np.sqrt(N).astype(int)\n",
    "    pxy,_,_ = np.histogram2d(x, y, bins)\n",
    "    pxx = np.tile(pxy.sum(1), (bins, 1)).T\n",
    "    pyy = np.tile(pxy.sum(0), (bins, 1))\n",
    "    mask = pxy!=0\n",
    "    mi_val = np.sum(pxy[mask]*np.log(pxy[mask]/pxx[mask]/pyy[mask]))\n",
    "    return mi_val / N + np.log(N)\n",
    "\n",
    "def dmi(x, y, delay):\n",
    "    if delay == 0:\n",
    "        return mi(x, y)\n",
    "    elif delay < 0:\n",
    "        return mi(x[-delay:],y[:delay])\n",
    "    elif delay > 0:\n",
    "        return mi(x[:-delay],y[delay:])\n",
    "\n",
    "def tdmi(x, y, time_range):\n",
    "    return np.array([dmi(x,y,delay) for delay in time_range])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "#### [Python code for mutual information estimation](https://github.com/NeoNeuron/mutual_information)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# install minfo, a Cython-based mutual information estimator;\n",
    "!pip install minfo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# import mutual information estimator\n",
    "from minfo.mi_float import tdmi_omp\n",
    "# wrap funciotns estimating time-delayed mutual information\n",
    "def TDMI(x:np.ndarray, y:np.ndarray, time_range:np.ndarray) -> np.ndarray:\n",
    "\n",
    "    assert 0 in delay\n",
    "    assert time_range.dtype == int\n",
    "\n",
    "    tdmi_val = np.zeros_like(time_range, dtype=float)\n",
    "    pos_end = np.max(time_range).astype(int)\n",
    "    neg_end = -np.min(time_range).astype(int)\n",
    "    if neg_end == 0:\n",
    "        tdmi_val = tdmi_omp(x, y, pos_end+1)\n",
    "    elif pos_end == 0:\n",
    "        tdmi_val = np.flip(tdmi_omp(y, x, neg_end+1))\n",
    "    else:\n",
    "        tdmi_val[neg_end:] = tdmi_omp(x, y, pos_end+1)\n",
    "        tdmi_val[:neg_end+1] = np.flip(tdmi_omp(y, x, neg_end+1))\n",
    "\n",
    "    return tdmi_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# test TDMI()\n",
    "test_x = np.random.randn(100000)\n",
    "test_y = np.random.randn(100000)\n",
    "test_y[1:] += test_x[:-1]*0.5\n",
    "delay = np.arange(21)-10\n",
    "plt.plot(delay, TDMI(test_x, test_y, delay))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "---\n",
    "# Section 2: Causal inference in linear and nonlinear systems.\n",
    "\n",
    "Now use these two methods you code above to play with few artifitially generated time series."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Case A: Unidirection misinferred as no interaction\n",
    "\\begin{array}{l}x_{t}=\\varepsilon_{t} \\\\ y_{t}=-0.1 x_{t-1}+\\eta_{t}\\end{array}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "def system_A(length):\n",
    "    x = np.random.randn(length)\n",
    "    y = np.random.randn(length)\n",
    "    y[1:] += -0.1*x[:-1]\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# GC version\n",
    "length = int(1e7)\n",
    "order = 10\n",
    "X, Y = system_A(length)\n",
    "print(f'GC Y->X     : {GC(X, Y, order):.3e}')\n",
    "print(f'GC X->Y     : {GC(Y, X, order):.3e}')\n",
    "print(f'GC threshold: {GC_SI(0.001, order, length):.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# TDMI version\n",
    "delay = np.arange(21)-10\n",
    "plt.plot(delay, TDMI(X,Y,delay))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Nonlinear transformation\n",
    "$$\n",
    "\\tilde{x}_t = x_t^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# GC version\n",
    "X_tilde = X**2\n",
    "print(f'GC Y->X     : {GC(X_tilde, Y, order):.3e}')\n",
    "print(f'GC X->Y     : {GC(Y, X_tilde, order):.3e}')\n",
    "print(f'GC threshold: {GC_SI(0.001, order, length):.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# TDMI version\n",
    "delay = np.arange(21)-10\n",
    "plt.plot(delay, TDMI(X_tilde,Y,delay))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Case B: Unidirection misinferred as bidirection\n",
    "\\begin{array}{l}x_{t}=-0.3 x_{t-1}+\\varepsilon_{t} \\\\ y_{t}=0.3 y_{t-1}-0.9 x_{t-1}+\\eta_{t}\\end{array}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "def system_B(length):\n",
    "    x = np.random.randn(length)\n",
    "    y = np.random.randn(length)\n",
    "    for i in range(1, length):\n",
    "        x[i] += -0.3*x[i-1]\n",
    "        y[i] += 0.3*y[i-1] - 0.9*x[i-1]\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# GC version\n",
    "length = int(1e7)\n",
    "order = 10\n",
    "X, Y = system_B(length)\n",
    "print(f'GC Y->X     : {GC(X, Y, order):.3e}')\n",
    "print(f'GC X->Y     : {GC(Y, X, order):.3e}')\n",
    "print(f'GC threshold: {GC_SI(0.001, order, length):.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# TDMI version\n",
    "delay = np.arange(21)-10\n",
    "plt.plot(delay, TDMI(X,Y,delay))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Nonlinear transformation\n",
    "$$\n",
    "\\tilde{x}_t = \\left[\\frac{x_t+|x_t|}{2}\\right]^5\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# GC version\n",
    "X_tilde = ((X+np.abs(X))/2)**5\n",
    "print(f'GC Y->X     : {GC(X_tilde, Y, order):.3e}')\n",
    "print(f'GC X->Y     : {GC(Y, X_tilde, order):.3e}')\n",
    "print(f'GC threshold: {GC_SI(0.001, order, length):.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# TDMI version\n",
    "delay = np.arange(21)-10\n",
    "plt.plot(delay, TDMI(X_tilde,Y,delay))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Case C: Unidirection misinferred as reversed unidirection\n",
    "\\begin{aligned}\n",
    "&x_{t}=-\\sum_{k=1}^{8} c_{k} x_{t-k}+\\varepsilon_{t} \\\\\n",
    "&y_{t}=-\\sum_{k=1}^{8} c_{k} y_{t-k}+100 \\sum_{k=1}^{9} c_{k-1} x_{t-k}+\\eta_{t}\n",
    "\\end{aligned}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "def system_C(length):\n",
    "    x = np.random.randn(length) * 2.27e-3\n",
    "    y = np.random.randn(length)\n",
    "    c_coef = np.array([-5.1777,12.6133,-18.6167,18.1332,-11.9147,5.1664,-1.3573,0.1678])\n",
    "\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# GC version\n",
    "length = int(1e7)\n",
    "order = 10\n",
    "X, Y = system_C(length)\n",
    "print(f'GC Y->X     : {GC(X, Y, order):.3e}')\n",
    "print(f'GC X->Y     : {GC(Y, X, order):.3e}')\n",
    "print(f'GC threshold: {GC_SI(0.001, order, length):.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# TDMI version\n",
    "delay = np.arange(21)-10\n",
    "plt.plot(delay, TDMI(X,Y,delay))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Nonlinear transformation\n",
    "$$\n",
    "\\tilde{x}_t = x_t^5\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# GC version\n",
    "X_tilde = X**5\n",
    "print(f'GC Y->X     : {GC(X_tilde, Y, order):.3e}')\n",
    "print(f'GC X->Y     : {GC(Y, X_tilde, order):.3e}')\n",
    "print(f'GC threshold: {GC_SI(0.001, order, length):.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# TDMI version\n",
    "delay = np.arange(21)-10\n",
    "plt.plot(delay, TDMI(X_tilde,Y,delay))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Case D: Bidirection misinferred as no interaction\n",
    "\\begin{array}{l}x_{t}=-0.1 y_{t-1}+\\varepsilon_{t} \\\\ y_{t}=-0.1 x_{t-1}+\\eta_{t}\\end{array}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "def system_D(length):\n",
    "    x = np.random.randn(length)\n",
    "    y = np.random.randn(length)\n",
    "    for i in range(1, length):\n",
    "        x[i] += -0.1*y[i-1]\n",
    "        y[i] += -0.1*x[i-1]\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# GC version\n",
    "length = int(1e7)\n",
    "order = 10\n",
    "X, Y = system_D(length)\n",
    "print(f'GC Y->X     : {GC(X, Y, order):.3e}')\n",
    "print(f'GC X->Y     : {GC(Y, X, order):.3e}')\n",
    "print(f'GC threshold: {GC_SI(0.001, order, length):.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# TDMI version\n",
    "delay = np.arange(21)-10\n",
    "plt.plot(delay, TDMI(X,Y,delay))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Nonlinear transformation\n",
    "$$\n",
    "\\tilde{x}_t = x_t^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# GC version\n",
    "X_tilde = X**2\n",
    "print(f'GC Y->X     : {GC(X_tilde, Y, order):.3e}')\n",
    "print(f'GC X->Y     : {GC(Y, X_tilde, order):.3e}')\n",
    "print(f'GC threshold: {GC_SI(0.001, order, length):.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# TDMI version\n",
    "delay = np.arange(21)-10\n",
    "plt.plot(delay, TDMI(X_tilde,Y,delay))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "## Case E: Bidirection misinferred as unidirection\n",
    "\\begin{aligned}\n",
    "x_{t}=-\\sum_{k=1}^{8} c_{k} x_{t-k}+0.5 \\sum_{k=1}^{9} c_{k-1} y_{t-k}+\\varepsilon_{t} \\\\\n",
    "y_{t}=-\\sum_{k=1}^{8} c_{k} y_{t-k}+0.5 \\sum_{k=1}^{9} c_{k-1} x_{t-k}+\\eta_{t}\n",
    "\\end{aligned}\n",
    "$$\n",
    "x_{t}=-\\sum_{k=1}^{8} c_{k} x_{t-k}+0.5 \\sum_{k=1}^{9} c_{k-1} y_{t-k}+\\varepsilon_{t} \\\\\n",
    "y_{t}=-\\sum_{k=1}^{8} c_{k} y_{t-k}+0.5 \\sum_{k=1}^{9} c_{k-1} x_{t-k}+\\eta_{t}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "def system_E(length):\n",
    "    x = np.random.randn(length)\n",
    "    y = np.random.randn(length)\n",
    "    for i in range(1, length):\n",
    "        # Add your code here\n",
    "        pass\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# GC version\n",
    "# length = int(1e2)\n",
    "# order = 10\n",
    "# X, Y = system_E(length)\n",
    "# print(f'GC Y->X     : {GC(X, Y, order):.3e}')\n",
    "# print(f'GC X->Y     : {GC(Y, X, order):.3e}')\n",
    "# print(f'GC threshold: {GC_SI(0.001, order, length):.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# TDMI version\n",
    "# delay = np.arange(21)-10\n",
    "# plt.plot(delay, TDMI(X,Y,delay))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Nonlinear transformation\n",
    "$$\n",
    "\\tilde{x}_t = \\tanh(10x_t)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# GC version\n",
    "# X_tilde = np.tanh(10*X)\n",
    "# print(f'GC Y->X     : {GC(X_tilde, Y, order):.3e}')\n",
    "# print(f'GC X->Y     : {GC(Y, X_tilde, order):.3e}')\n",
    "# print(f'GC threshold: {GC_SI(0.001, order, length):.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# TDMI version\n",
    "# delay = np.arange(21)-10\n",
    "# plt.plot(delay, TDMI(X_tilde,Y,delay))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "# Summary Table\n",
    "\n",
    "![](./static/summary_table.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "---\n",
    "# Section 3: Interactions between neuronal signals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "from scipy.io import loadmat\n",
    "try:\n",
    "    data = loadmat('lab02-xk128-103005006-f.mat')\n",
    "except FileNotFoundError:\n",
    "    data = loadmat('./tutorials/W6_NetworkCausality/lab02-xk128-103005006-f.mat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "spk_train = data['scsig053ats'][0,0][-1]\n",
    "lfp = data['AD27_ad_000'][0,0][-1]\n",
    "\n",
    "t_start = 170   # unit ms\n",
    "t_end   = 270   # unit ms\n",
    "\n",
    "dt = 0.001\n",
    "lfp_theta = lfp[int(t_start/dt):int(t_end/dt)].flatten()\n",
    "spike = spk_train[(spk_train>=t_start)&(spk_train<t_end)] - t_start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "def spike2bool(spike_train, length, dt):\n",
    "    '''\n",
    "    Transform spike train data to binary time series.\n",
    "\n",
    "    Args:\n",
    "    spike_train : original spike train\n",
    "    length      : length of desired binary series\n",
    "    dt          : size of time step\n",
    "\n",
    "    Return:\n",
    "    binary time series of spiking time series\n",
    "    '''\n",
    "    spike01 = np.zeros(length, dtype=int)\n",
    "    spike01[np.floor(spike_train/dt).astype(int)] = 1\n",
    "    return spike01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "spike01 = spike2bool(spike, len(lfp_theta), dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "delay = np.arange(1001)-500\n",
    "plt.plot(delay, tdmi(spike01, lfp_theta, delay), 'navy')\n",
    "plt.axvline(0, color='orange')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# GC analysis\n",
    "order = 30\n",
    "print(f'GC lfp->spike : {GC(spike01, lfp_theta, order):.3e}')\n",
    "print(f'GC spike->lfp : {GC(lfp_theta, spike01, order):.3e}')\n",
    "print(f'GC threshold  : {GC_SI(0.001, order, len(lfp_theta)):.3e}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "# End of the tutorial\n",
    "---\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W6_Tutorial2",
   "toc_visible": true
  },
  "interpreter": {
   "hash": "9516f62da91337f10c2adbe814d9c63a4b08f8271333386358218606edb781e3"
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3.7.11 64-bit ('pw3': conda)",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
